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ABSTRACT Given two sets of functional data having a common underly¬ 
ing mean function but different degrees of distortion in time measurements, 
we provide a method of estimating the time transformation necessary to align 
(or ‘register’) them. We prove that the proposed method is consistent under 
fairly general conditions. Simulation results show superiority of the perfor¬ 
mance of the proposed method over two existing methods. The proposed 
method is illustrated through the analysis of three paleoclimatic data sets. 
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1. Introduction 

Gonsider functional data arising from observations recorded at a sequence 
of time points. The task of aligning multiple but similar sets of functional 
data by possibly nonlinear adjustment to their time scales is often referred 
to as ‘registration’. The problem of registration is also important in image 
and video processing, where multiple dimensions are involved. In the one 
dimensional case, the need for registration has been felt in broadly two types 
of problems. For longitudinal growth data, often viewed as a common pattern 
expressed differently through different individuals with their diverse scales of 
evolution, the need for registration arises from the quest for the common 
pattern. In this type of problems, the number of individuals is generally 
much more than the number of observations per individual. On the other 
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hand, for paleoclimatic data on historic movement of climate variables, the 
need for registration arises from the notion that the recorded ‘time’ contains 
estimation error that might come in the way of collation of information from 
multiple sets of data. In these applications, the number of observations per 
data set is much more than the number of data sets to be time-aligned or 
registered. There may even be only two data sets for alignment. In this 
paper, we address the second type of problems, and focus on registration of 
one data set with respect to another. While registration may be needed for 
the purpose of comparing/correlating one variable with another, we consider 
it only in the context of pooling two data sets on a common variable for the 
ultimate purpose of better description of its movement over time. 

As a motivating example, consider the atmospheric concentration of car¬ 
bon dioxide for the past 400,000 years constructed from two ice cores ex¬ 
tracted from two different locations in Antarctica, plotted in Figu re [H One o 
these data series was collected from the ice sheet over lak e Vos tok flPetit et al. 


19991) . and the other from EPICA dome flLiithi et al. 


20081). The peaks in 


the EPICA dome data (many of which are sharp) precede those in the lake 
Vostok data during the initial and later parts of the time span, while the 
reverse happens in the middle part. It suggests the possible existence of a 
non-linear relationship between the time values of the two data sets. 

Figure [U here. 


Sharp peaks and valleys present in functional data sets may be viewed 
as either a help or a hindrance for registration. One may hold the view 
that paying too much attention to these details may cause distortion, as ob¬ 
servation times in either of the data sets may miss the actual peaks of the 
underlying continuous time phenomenon. In that case, the method of regis¬ 
tration should not utilize these features as such. Several methods of this kind 


are available: contin uous monotone reg i stration flRamsav and Lil. Il998h . dy- 


gression flKneip et al 


namic ti me warning fWang and Gasserl . 119971.1199911 . registration by local re- 

ration through para- 


20001) , maximum like 


metr ic modeling of the time tran sformation (|R,Onnl. 


ihooc 


regis 


2001 


20051) . self-modelling warping flGervini and Gasserl . 1200411 . shape invariant 


Gervini and Gasserl. 
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model based registration flRrnmback and Lindstroml. 


2004), fnnctional con¬ 


vex synchronization model based registration flLin and Miillerl. l2004l) . pair 


wise curve synchronization flTang and Mill 


erl. 


1995 


20081). and functional princip al 


Kneip and Ramsavl . 120081) . 


component based registration flSilvermanl . 

One may also hold the view that sharp peaks and valleys in the data are 
characteristics of the underlying common function, and in that case it may 
be unwise to ignore the information contained in them. Making use of such 
information may yield better registration. Smoothing based methods may be 
unable to take full advantage of this information, as smoothing blurs these 
features. 

It is possible to use only the sharp peaks and valleys in the data series, 
often identihed as ‘landmarks’, for the purpose of registration. A simplistic 
method would be to identify specific landmarks of one data set that cor¬ 
respond to the landmarks of the other data set, and use a piecewise linear 
time transformation th at permits the alignment of the two sets of landmarks. 


Ramsay and Lil f 19981) suggested that such a crude registration may be fur¬ 
ther rehned through their method. The other methods mentioned above 
may also be used for rehnement. Manual identihcation of some matching 
landmarks in the two data sets (after preli minary smoothing) have been 
used as input to the m ethods proposed bv iKnein and Gasser! fll992l) and 


Kneip and Engell (119951) . The method of iKneip and Gasserl (119921) leads to 


registration under the constraint that the identihed landmarks should match. 
The identihcation of matching landmarks is generally done manually. This 
is a disadvantag e of landmark-assisted registration. 

Bigotl (120061) proposed an automated landmark-based method, which pro¬ 
vides registration through the following steps: (a) identihcation of signih- 
cant landmarks in each of the two data sets, (b) establishing possible corre¬ 
spondence between these landmarks, and (c) estimation of the time trans¬ 
formation function through a standard method of nonparametric smooth¬ 
ing/regression, by using the pairs of matched times of landmarks as input. 
The performance of this multiple-step method may be limited by the non-use 
of information other than landmarks and possible accumulation of estimation 
errors at diherent steps. 
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(which is natural in longitudinal growth data) for estimation of the popula¬ 
tion characteristics, and are not usable for applications where only a pair of 
data sets need to be registered. 


With the objective of registering one functional data set with respect to 
another, we model the observations of the two data sets as 

yi{t) = m{t) + ei{t) 
y^{t) = rn{g{t)) + e 2 {t) 

where m is the underlying mean function, g is the requisite time transfor¬ 
mation, and Ci and 62 are random error terms. The task of registration is 
essentially that of estimating g that aligns the two sets of data. Note that the 
above model would not be suitable for data sets representing observations on 
different variables (e.g., paleoclimatic data sets on temperature and carbon 
dioxide). On the other hand, any generalization of the above model for ap¬ 
plicability to different variables would involve additional estimation, which 
is redundant when the observed variables are in fact identical. Thus, such 
a generalization can only be achieved at the cost of inefficient utilization of 
the available information. 

In Section 2, we propose a new estimator of the time transformation 
function g^ by maximizing a ‘measure of alignment’ of two functional data 
sets. The maximization is done over an appropriate class of transformations. 
This measure of alignment is designed to capture the information contained 
in the entire set of data, including locations of sharp variation. We show 
that it possesses some desired characteristics. The method is automated as 
it does not require manual identihcation of landmarks. Identihability of g 
with respect to the chosen model, under appropriate conditions, is proved in 
Section 3. 

Many of the existing methods of estimation of the time t ransformation 
have been proposed without establishing their consistency. 


R.onnl ( 1200 ih 
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and 


Gervini and Gasser! (120051) have proved the consistency of their methods 


when t he number of subj e cts go es to inhnity. Following 


fll995fi . IWang and Gasser 


fll999h . and 


Gervini and Gasser 


Knein and Engei 


fl2004h . we have 


established the consistency of our estimator as the numbers of observations 
in the two data sets go to inhnity. This result, reported in Section 3, holds 
when the time transformation is chosen from a class of functions satisfying 
some general conditions. 

Results of a simulation study, to demonstrate the performance of the 
estimator chosen from a particular class, are reported in Section 4. The 
method is illustrated in Section 5 through the analysis of seyeral real data 
sets. Some concluding remarks are provided in Section 6 . 


2. Model and Methodology 

Let and {(si, R/J,..., (s^^, be two sets of 

functional data, arising from the model 

Yti = m{ti)+ei{ti) i = l,...,ni 

= "i(^o(si))+ e 2 (sj), j = l,...,n 2 ( 1 ) 


where m is an underlying location function that is continuous, and qq is 
an unknown time transformation function, which is continuous and strictly 
increasing. The terms Ci and 62 represent additive random measurement 
errors, which have mean zero. 

Now, let us dehne, for any given continuous and strictly increasing trans¬ 
formation function g, the functional 



where n = rii -|- 77 - 2 , Ki and K2 are kernel functions that are probability 
densities, and hi and /i 2 are the corresponding bandwidths. The above func- 

( Yj, _y' \ 

i = 1 ,... ,ni and j = 1,..., n2. The weights depend on g. Note that when 
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g = go, for every pair of i and j such that g{sj) is close to tj, the continuity of 
m ensures that ht. — is expected to be small. Therefore, for g = go, large 

( Vt ■ —Y' \ 

j are expected to occur together with large values of 
their weights. This may not be the case when g ^ go- In Section 3, we show 
that under some general conditions, the probability limit of Ln{g) attains its 
maximum value if and only ii g = go- 

Thus, Ln{g) may be interpreted as a measure of alignment. 

Let us now examine the roles of the bandwidth parameters hi and h 2 in 
the above measure. A small value of hi makes the weight for a given i and j 
nearly equal to zero, unless U is very close to g{sj). Thus, only a few weights 
can be substantial. When hi is large, weights can be substantial for more 
combinations of i and j. Thus, hi controls the effective number of weights 
in the weighted sum in ([2]). On the other hand, /i 2 controls the penalty 
for discrepancies between and Y^.. A very large value of ^2 might make 
Ln{g) insensitive to changes in g, as there would not be enough penalty for 
mismatch between Yt^ and Y^.. A very small value of h 2 would make Ln{g) 

( Yt . —Y' \ 

j would be nearly zero for most of the combinations 

of i and j. 

We dehne the proposed estimator of the function go as 

gn = aigmi^Lnig), (3) 

geG 

where Ln{g) is as dehned in (|2]) and G is a suitable class of continuous and 
strictly increasing functions that includes the true transformation function go- 
As for choices of the bandwidths, one can select hi as a fraction of the 
range of time in either sample and /i 2 as a fraction of the combined range of 
the observed variable. Some guidelines are given in Section 4. 

Corresponding peaks in two sets of functional data are sometimes iden- 
tihed manually as matching landmarks. In our case, the objective function 
Ln{g) automatically rewards candidate transformation functions that map 
peaks of one data set into the corresponding peaks of the other. On the 
other hand, if a peak is missing from one of the data sets, then it does not 
penalize the ‘correct’ transformation any more than a similar alternative can- 
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didate (i.e., a marginally different transformation). Therefore, the estimator 
Qn should be able to utilize the landmarks automatically for registration. 


3. Consistency 

Let the errors {ei(tj); i = 1,... rii} and { 62 ( 5 ^); j = 1,... n 2 } and the time 
points {ti] i = l,...ni} and {sj] j = l,...n 2 } be mutually independent 
sets of samples from the probability density functions fe^, /ej, fi, and /2 
having supports over [—cxo, cxo], [—cxo, cxo], [a, b] (for 0 < a < 6), and [c, d] (for 
0 < c < d) respectively. 

We assume that m in model m is a continuous function dehned over 
[0, cxo). We also presume that it is not flat or it does not fluctuate too much. 
In order to ensure this, we formally stipulate that for any given interval [p, q ], 
the inverse image (with respect to m) of any point on m([p, g]) has a hnite 
intersection with [p, g], i.e., for every y G [minig[p^q] m(t), maXig[p^q] m(t)], the 
set {t : m{t) = y, t E [p, g]} has a hnite number of elements. 

Estimation of go ia model ([I]) makes sense only if go{[c, d\) has a substan¬ 
tial overlap with [a, b]. We, therefore, assume that go{[c, d]) fl [a, b] includes a 
non-empty open interval. We also need to ensure that there is no ambiguity 
about go in model ([T]). Let Go be the class of all strictly increasing and con¬ 
tinuous functions g dehned over [c, d] such that the set Sg = g~^{[a, 6]) fl [c, d] 
contains a non-empty open interval and that g agrees with go at least at one 
point in Sg^ H Sg. By construction, go ^ Go. Our hrst result establishes the 
identihability of go (within Go) with respect to model ([1]). 


Theorem 3.1 Let m, go, and Go be as described above. If g E Go is such 
that m(g(s)) = m(go(s)) for all s E SgO then Sg = Sg^ and g(s) = go(s) 
for all s E Sg^. 


The consistency of g^ n eeds to be establis 
the data sets go to inhnity flKneip and Engei 


red as the sample size in both 


1995 


Wang and Gasserl. 


1999 


Gervini and Gasserl. 120041) . As a hrst step, we establish the point-wise con¬ 


vergence of the functionals Ln on Gq after making the following assumptions. 


Assumption 1. The densities f^^, /i, and /2 are continuous and 
bounded; and f^.^ symmetric about zero and are strictly unimodal 
at zero; fi and /2 are positive over the interior of their supports. 
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Assumption 2. The kernels Ki and K 2 are continuous and bounded prob¬ 
ability density functions defined over the real line. 

This condition is satished by all the popular kernels viz. Uniform, Triangular, 
Epanechnikov, Biweight, Gaussian, and so on. 

Assumption 3. The sample sizes rii and 77-2 are such that ni/n —)■ ^ for 
some ^ G (0,1), as n —)■ 00 . 

Assumption 4. The bandwidths are such that hj —)■ 0 and njhj —)■ cxo as 
n — 00 , 7 = 1, 2. 


Theorem 3.2 Let m, go, and Gq be as described at the beginning of this 
section. Then, under Assumptions I- 4 , for any function g G Go, as n ^ 00 , 

p 

Ln{g) L{g), where 

T( \ f^^9{.y))f2{,y)feAv-^i.9i.y)) + 'm{go{y)))fe 2 {v)dydv 

- FM9{y))My)dv ' ''' 

We now show that this limiting functional is maximized only by the correct 
transformation function. 


Theorem 3.3 Suppose m, go, and Gq are as described at the beginning of 
this section. Then, under Assumption 1, 

(a) L{g) < L{go) for all g G Go, 

(b) If L{g) = L{go) for some g G Go, then g = go over Sg^. 

The next step is to establish the uniform convergence of L„, for which we 
need a stronger condition on Go that enforces compactness. Let us dehne a 
metric on Go viz., d{gi, g 2 ) = sup^gj^^^j \gi{x)-g 2 ix)\ = || 5 'i- 5 ' 2 ||; 91,92 e Go. 

Assumption 5. The class G in ([3]) is a compact subset of Go in the metric 
space (Go, d) and it includes go- 

An example of G that satishes Assumption 5 is the subset of functions g of 
Go with bounded slope. 
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Assumption 2A. The kernels Ki and K 2 are bounded away from zero on 
a given closed interval and have bounded hrst order derivatives. 

The Gaussian kernel satishes the above condition. 

Theorem 3.4 Let m, go, and Go be as deseribed at the beginning of this 
section. Then under the Assumptions 1-5, and 2A, as n ^ 00 , 

p 

sup \Ln{g) - L{g)\ 0. 

g€G 

We now establish that the sequence of maximum values of the functionals 
Ln converges to the value of L at its maximizer, go. 

Theorem 3.5 Let m, go, and Gq be as described at the beginning of this 

p 

section. Then under the Assumptions 1-5, and 2A, Ln{gn) —^ L{go) as 
n ^ 00 . 

Finally we establish the consistency of our estimator. 

Theorem 3.6 Let m, go, and Gq be as described at the beginning of this 

p 

section. Then under the Assumptions 1-5, and 2A, gn ^ go as n ^ 00 . 

4. Simulation of performance 
4.1.Methods compared 

Even though we have proposed a general class of estimators in Section 2 
and established their consistency in Section 3, we need to focus on a specihc 
member of that class in order to simulate the performance. We chose G as 
the vector space generated by linear B-spline basis functions with equidistant 
knot points over [a,b]. The method of Steepest Ascent was used to maxi¬ 
mize (j2]). We opted for standard normal densities for Ki and K 2 . We chose 
hi as 5 per cent of the range of time in either data set, and h 2 as 10 per cent 
of the range of combined T-values of both the data sets. In order to provide 
an initial iterate to steepest ascent, we programmatically mapped signihcant 
peaks and valleys of the two data sets. The optimization algorithm, however, 
produced reasonably good results (not reported here) even when the identity 
map was provided as initial iterate. 
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We carried out simulations to compare the performance of the above 
implementation of the proposed method with two other methods. 

1. The hrst method was continuous monotone registration. We used its 
implementation in the R-function register.fd in the package fda. We 
retained default values of the order of the polynomial splin es (4) and the 
roughness penalty (2). Following Ramsav and Li fll998 h we selected 
the value of the smoothing parameter A as 10“^. 

2. The second method is self-modelling registration. We used the matlab 
codes provided by Daniel Gervini on his web site (https: //pantherf ile. 
We chose default values for the number of random starts (2 0) and the 
order of splines (3). Following iGervini and Gasser! (120041) other pa¬ 
rameters of the method viz. the numbers of components (g) and the 
number of basis functions (p) were selected by using the cross-validation 
algorithm. 

These methods were selected for comparison mostly on account of availability 
of codes and applicability to the data at hand. The automated method of 
Bigotl (120061) was not chosen for comparison, as his code can only handle 
sample sizes that are powers of 2. 


4.2. Simulation Design 

The choice of data for simulation was motivated by the paleoclimatic data 
set that exhibits sharp changes. We chose lake Vostok data on atmospheric 
concentration of carbon dioxide (with 283 data points, range of time-values 
/ = [2.3,414.1] and s. d. of R-values as s = 28.7) described in Section 1, as 
the base for our simulation exercises. We conducted Monte Garlo simulations 
in four different scenarios under model ([T]) (with rii = n 2 = 250) as described 
below: 


Scenario 1. Here we randomly selected a sub-set D of 250 data points from 
the base set. The function m was obtained by linear interpolation from 
those selected points. In every simulation run, we constructed the pair 
of data sets as follows. 
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(a) Time values of both the data sets were kept same as those of D. 

(b) The time transformation qq in ([1]), applied on the time values 
of the second data set, was chosen as the linear spline with six 
eqnidistant knots over the range of time-valnes of D (see graph in 
Figure ??(i) in Supplemental Materials), i.e.. 


g^it) = -0.379 + 1.05f- 0.209(t-85.7)+ 

+ 0.409(t - 167.8)+ - 0.609(f-249.9) + 

+ 0.809(t-332.0)+ , for t e [3.6,414.1] (5) 


where 



u if M > 0, 
0 if M < 0. 


(c) Additive errors were generated afresh from the normal distribution 
(mean=0, s.d.=0.05s) for both the data sets separately. 

Scenario 2. Here, m was the fnnction obtained by linearly interpolating 
between the data points of the base data. The 250 time points as well 
as the additive errors for each data set were generated afresh for each 
simnlation run. In this case, the distributions fi and /2 described in 
Section 3 were uniform over I (see the beginning of this snbsection), 
while the distribntions were normal (mean=0, s.d.=0.05s). The 

function go was as in ([5]). 

Scenario 3. This set-up was the same as in Scenario 1, except that go was 
chosen as the identity fnnction plus a periodic function viz.. 


go{t) = t + 0.05t sin(37rt/&) for t E [a, b]. 


( 6 ) 


The graph of this function is shown in Figure ??(ii) in Supplemental 
Materials. 

Scenario 4. This set-np was the same as in scenario 2, except that go was 


as m 
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Scenario 1 offers an opportunity to assess the performance of the methods 
when G includes go- We chose 21 equidistant knots in our search space of 
B-splines, so that go with six equidistant knots becomes a special case. In 
Scenarios 2, 3, and 4, the chosen number of knots was 20, which means that 
go is not included in the search space. For continuous monotone registration, 
the number of knots chosen for each scenario was the same as that of the 
proposed method. 

The matlab programs for self-modelling registration require values of the 
functions to be registered at a common set of time points, which is met 
only in scenario 1 and 3. Therefore, simulation results for this method are 
reported only for these two scenarios. 

The performance of the estimators of go were studied in terms of (a) point- 
wise bias, (b) point-wise standard deviation, (c) point-wise mean squared 
error (MSE), and (d) average of the integrated mean square error (IMSE) 
normalized by the squared norm of the true function, dehned for each simu¬ 
lation run as 

SiaKtidt 

where S is the number of independent runs of the simulation and gj is the 
estimate of go at the jth run. We used Simpson’s rule to evaluate these 
dehnite integrals. 

4.3. Results 

The results of the simulations based on 1000 independent runs for each of 
the scenarios, are shown in Figure [2l In scenarios 2 and 4, where time- 
points change from one run to another, the R-function register.fd did not 
produce results for 240 simulation runs. Therefore, results for this estimator 
corresponding to these scenarios are based on 760 runs. 

Figure [2] here. 

Table [U here. 

It is observed that, the proposed estimator has, in general, smaller bias, 
standard deviation, and MSE as compared to continuous monotone registra¬ 
tion, and self-modelling registration. However, towards the right end of the 
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time scale in Scenarios 3 and 4, the proposed method exhibits both higher 
standard deviation, and MSE as seen in Figure [2l In this context, it may be 
mentioned that Kernel based methods are known to perform poorly near the 
fringes of the data. It is observed in Table [T] that the proposed method had 
uniformly smaller average normalised IMSE for all the scenarios. 


5. Analysis of ice core data 

Here we considered paleoclimatic d ata on the a 


carbon dioxide, and (ii) methane fjPetit et al. 


mospheric concentration of fi) 


I 999 I: Ihoulergue et all. 120081) 


as determined from air-bubbles trapped in ice cores collected over Lake Vos- 
tok and at EPICA Dome of Antarctica, and fiii) average annual temperature 


deviations fjPetit et al. 


1999 


Masson-Delmottel . 120071) . which were recon¬ 


structed from deuterium contents at various depths of ice cores obtained at 
these two sites. Table [2] gives some descriptive statistics of the data sets. 


Table |2] here. 


We chose to align the data set from EPICA dome with that from lake 
Vostok using two registration methods viz. the proposed method, and contin¬ 
uous monotone registration. We could not apply the method of self-modelling 
registration for the above data, since its matlab implementation requires the 
nominal observation times in the two data sets to coincide. 

The number of knots used for both the methods compared here were 30, 
30, and 20 for carbon dioxide, methane, and temperature deviation data 
respectively, while remaining parameters were chosen as in the previous sec¬ 
tion. 

Alignments of the carbon dioxide, temperature, and methane data sets 
are shown in Figures [21IH and |2] respectively. The R-function register.fd 
could not produce output for the methane data. The estimates of the function 
go for the three pairs of data sets, produced by the methods are plotted in 
Figure ?? of Supplemental Materials. 

Figure |3] here. 


Figure m here. 
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Figure O here. 

Table [3] here. 

The average squared distance between the interpolated curves before and 
after registration by the two methods, computed over a uniform grid of size 
1000 over the common time range, are reported in Table [3l It is observed 
that the proposed method produced shorter distance between the registered 
curves as compared to continuous monotone registration. 

5. Discussion 

In this paper, we have proposed a new method of registration of one func¬ 
tional data set with another, by optimizing an empirical measure of align¬ 
ment. If there are sharp variations in the data, the proposed method is able 
to utilize them, without requiring prior identihcation of landmarks. Since 
the method does not use any pre-smoothing, it does not suffer from any loss 
of information that might occur due to smoothing. On the other hand, the 
measure of alignment (2) ensures that the proposed method makes use of the 
main strength of smoothing, namely pooling of information from neighbour¬ 
ing observations. 

The present implementation of the method, in the form of an R code, is 
available from the authors on request. This implementation permits registra¬ 
tion of data sets with possibly unequal, irregularly spaced and large number 
of samples. This implementation is based on some specihc choices, e.g., use 
of the class of B-splines with uniformly spaced knot points as candidate time 
transformation functions, and steepest ascent for optimization. However, 
none of these choices is necessary in the general set-up used for proving the 
consistency of the proposed class of estimators. 

There are indeed some limitations of the proposed approach. The align¬ 
ment provided by the proposed method is likely to change if the data sets 
for registration are interchanged. When there are many sharp changes in the 
data, the iterative algorithm may run into spurious local maxima, particu¬ 
larly when the bandwidth parameter hi in (2) is small. Such problems may 
be mitigated by replacing steepest ascent with a probabilistic search algo¬ 
rithm such as simulated annealing. On the other hand, when the data sets 
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do not have sharp changes, the proposed method may not be very sensitive 
to alterations in the transformation function. The guideline on the choice 
of bandwidths (particularly /i 2 ) given here may not be applicable for those 
problems. The proposed method would be completely unsuitable for lon¬ 
gitudinal data with many individuals but relatively fewer observations per 
individual. 

When there are more than two data sets to be registered, the proposed 
method has to be used multiple times on pairs of data, possibly after iden¬ 
tifying one of the data sets as reference for registration. This reference data 
set may also be selected on a trial basis, and the candidate leading to the 
best overall alignment may be selected as reference data set at the end. 

The proposed method can be used as a tool for structural averaging 
(which amounts to estimating the function m in model (1)). Large sam¬ 
ple properties of the resulting estimator of m and its performance in relation 
to competing methods need to be studied in future. 

Registration of two sets of functional data on different variables (e.g., 
paleoclimatic data on temperature and carbon dioxide) is sometimes needed 
for the purpose of studying the relationship between them. The method pre¬ 
sented here is not readily applicable to this problem. However, an adaptation 
may be possible by inserting an unknown amplitude parameter in one of the 
two equations of model (1). This problem is another possible direction of 
future work in this area. 
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Figure SI. Qq functions for simulations 


Figure S2. Estimates of Qq for paleoclimatic data sets 


16 


Feature sensitive automated registration 


References 

Bigot, J. (2006). Landmark-based registration of curves via the continuous wavelet 
transform. J. Comput. Graph. Statist., 15, 542-564. 

Billingsley, P. (1985). Probability and Measure. John Wiley. 

Brumback, L. C. & Lindstrom, M. J. (2004). Self modeling with flexible, random 
time transformations. Biometrics, 60, 461-470. 

Gervini, D. &: Gasser, T. (2004). Self-modelling warping functions. J. R. Stat. 

Soc. Ser. B. Stat. MethodoL, 66, 959-971. 

Gervini, D. & Gasser, T. (2005). Nonparametric maximum likelihood estimation 
of the structural mean of a sample of curves. Biometrika, 92, 801-820. 

James, G. M. (2007). Gurve alignment by moments. Ann. Appl. Stat., 1, 480-501. 

Kneip, A. & Engel, J. (1995). Model estimation in nonlinear regression under 
shape invariance. Ann. Statist., 23, 551-570. 

Kneip, A. & Gasser, T. (1992). Statistical tools to analyze data representing a 
sample of curves. Ann. Statist., 20, 1266-1305. 

Kneip, A. &: Ramsay, J. O. (2008). Gombining registration and fitting for functional 
models. J. Amer. Statist. Assoc., 103, 1155-1165. 

Kneip, A., Li, X., MacGibbon, K. B., & Ramsay, J. O. (2000). Gurve registration 
by local regression. Canad. J. Statist., 28, 19-29. 

Liu, X. & Muller, H.-G. (2004). Functional convex averaging and synchronization 
for time-warped random curves. ,J. Amer. Statist. Assoe., 99, 687-699. 

Loulergue, L., Schilt, A., Spahni, R., Masson-Delmotte, V., Blunier, T., Lemieux, 

B., Barnola, J.-M., Raynaud, D., Stocker, T., , and Ghappellaz, J. (2008). 

Orbital and millennial-scale features of atmospheric ch4 over the past 800,000 
years. Nature, 453, 383-386. 

ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/antarctica/epica_doinec/edc-ch4-20 


Bhauniik et al. 


17 


Liithi, D., Floch, M. L., Bereiter, B., Blunier, T., Barnola, J. M., Siegenthaler, 

U., Raynaud, D., Jouzel, J., Fischer, H., Kawamura, K., &: Stocker, T. (2008). 

High-resolution carbon dioxide concentration record 650,000-800,000 years 
before present. Nature, 453, 379-382. 

ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/antarctica/epica_doniec/edc-co2-2008.txt 


Masson-Delmotte, V. (2007). Lsce/ipsl igbp pages/wdca contribution. (2007-091). 

ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/antarctica/epica_doniec/edc3deuttemp2007 

Petit, J. R., Jouzel, J., Raynaud, D., Barkov, N. I., Barnola, J. M., Basile, I., 

Bender, M., Chappellaz, J., Davis, M., Delaygue, G., Delmotte, M., Kotlyakov, 

V. M., Legrand, M., Lipenkov, V. Y., Lorius, C., Pepin, L., Ritz, C., Saltzman, 

E., &: Stievenard, M. (1999). Climate and atmospheric history of the past 
420,000 years from the vostok ice core, Antarctica. Nature, 399, 429-436. 

ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/antarctica/vostok/co2nat.txt,ftp://ftp. 

Ramsay, J. O. &: Li, X. (1998). Curve registration. J. R. Stat. Soc. Ser. B. Stat. 

MethodoL, 60, 351-363. 

Rpnn, B. B. (2001). Nonparametric maximum likelihood estimation for shifted 
curves. J. R. Stat. Soc. Ser. B. Stat. MethodoL, 63, 243-259. 

Silverman, B. W. (1995). Incorporating parametric effects into functional principal 
components analysis. J. R. Stat. Soc. Ser. B. Stat. MethodoL, 57, 673-689. 

Tang, R. Sz Muller, H.-G. (2008). Pairwise curve synchronization for functional 
data. Biometrika, 95, 875-889. 

Wang, K. &: Gasser, T. (1997). Alignment of curves by dynamic time warping. 

Ann. Statist, 25, 1251-1276. 

Wang, K. (fe Gasser, T. (1999). Synchronizing sample curves nonparametrically. 

Ann. Statist., 27, 439-460. 

Dibyendu Bhaumik, Department of Statistics and Information Management, Re¬ 
serve Bank of India, C9, 6th Floor, Bandra-Kurla Complex, Bandra (East), Mum¬ 
bai - 400051, India. 

E-mail: dbhaumik@rbi.org.in 



18 


Feature sensitive automated registration 


A Appendix: Proof of theoretical results 


Proof of Theorem 13.11 Let S' G Go be such that g{a) = go{a) for some a G 
^90 n Sg, m{g{s)) = m{go{s)) for all s G Sg^ fl Sg, and yet g{s) / go{s) for some 
s G SgQ n Sg. We can presume, without loss of generality, that g{s) < go{s). Let 
us assume, for now, a < s. Then the set 

{t: m{t) = m{go{s))}n[go{a),go{s)] 

has at least two elements, g{s) and go{s). Let < ••• < be the ordered 
elements of this set. Clearly, go{s) = Let g{s) = for some i < k. 

In order that the functions m{g{u)) — m{g{s)) and m{go{u)) — m{go{s)) coincide 
for all u G [a, s], these functions should have exactly the same number of zero 
crossings over this interval. However, from what we have already observed, the 
first function has exactly i zero-crossings, while the second function has exactly k 
zero-crossings, and i < k. Therefore, the two functions must differ somewhere on 
[a, s]. Similarly, if a > s, the set, 

{t: m{t) = m{g{s))}n[g{s),g{a)] 

has at least two elements viz., g{s) and go{s). If < • • • < be the ordered 
elements of the set, then g(s) = and go{s) = for some j > 1. Following 
similar arguments, the two functions above, which have been presumed to coincide 
for all u G [s,a], must differ somewhere on [s,®] as the first function has exactly I 
zero-crossings, while the second has exactly I — j + 1 . 

This contradicts the presumption that g{s) / go{s) for some s G Sg H Sg^. 
Thus, g = go over Sg n Sg^. 

The continuity of g and go, together with their equality over Sg fl Sg^, implies 
that Sg = SgQ. This completes the proof. □ 


Proof of Theorem 13.21 In accordance with Assumption 1, we denote by Mj- a 
common upper bound of the densities fi and / 2 . 

For a given time transformation g G Go, from m, we have 


Fnid) 


Nnjg) 

DnigY 


(7) 
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where 


Nn{g) = 


Dn{g) = 


T.T.kK. 


(ti - g{sj)\ 1 


nin2 ^ ^ hi' 

1=1 ]=i 


hi 


h2 


K2 


Y,, - Y' 


n\ n2 


1 1,, (u-aM 


nin2 ^ ^ hi 

1=1 j=i 


hi 


( 8 ) 

(9) 


We first establish that 

r*oo roo 


Nn{g)^ / fiig{y))f2{y)fei{v-m{g{y))+m{goiy)))fe2iv)dydv (10) 

Dn{g)^ [ fi{g{y))f2{y)dy. (ll) 

Jo 

The proof is then compl eted by using the continuous mapping theorem of conver¬ 
gence in probability fsee iBillingslevI (jl 985 l ll. 

Note that, from Q, we have 


EWnig)] = 


1 


nin2 


ni n2 

EE^ 


1 -g{sj)\ 1 


hi 


hi 


i=i j=i 

From the description of the model dU), we have 

'a; - g{y) 


■pK2 

h2 


'Y_^ 

h2 


/ OO roo poo poo -j 
-oo J—oo J ^ 


X ^K 2 

h2 


hi 

m{x) - m{go{y)) +u-v 

h2 


/i (a:) /2 (y) /ei («) /62 {v)dxdydudv. 


By making the transformations zi = ^ and Z2 = ^i9Y.y))+u v ^ 


h2 


have 


where 


/ OO poo poo poo 

/ / / Sn{zi,Z2,y,v)dzidz2dydv 

-oo^O J—oo J—oo 


Sn{zi,Z2,y,v) = I(^_g(^yyhi,oo){zi)Ki{zi) K2 {z 2) fi{g{y) + Zihi) 

^f2{y)fei{v - m{g{y) + zihi) + m{gQ{y)) + Z2h2)fe2{v)- 


As 5 is a positive and increasing function, any given real zi is contained in 
{—g{y)/hi,(x>) for sufficiently small hi. By using Assumptions 1 , 3 and 4 and 
the fact that m is a continuous function, for any fixed {zi, Z2,y,v), we have 

lim 5 ’n( 2 ;i,^ 2 , 2 /,i^) = Ki{zi)K2{z2) fi{g{y))f2{y) 

n—)-oo 

X fei{v-m{g{y))+m{go{y)))fe 2 iv). 


(12) 
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Note that, from Assumption 1 , we have 

0 < Sn{zi,Z 2 ,y,v) < MfKi{zi)K 2 {z 2 )f 2 {y)fe 2 {v)- (13) 

Assumption 2 ensures that the bounding function on the right hand side of (I13p is 
integrable. Then, by using the Dominated Convergence Theorem (DCT), we have 

/ oo roo 

/ fi{9iy))f2{y) 

-oo Jo 

X feiiv-m{g{y))+ m{go{y)))fe 2 {v)dydv. (14) 


From Lemma lA.li proved below, we have lim^^co Car[A„(o')] = 0. This estab¬ 
lishes ([To]) . 

We now turn to Dn{g). By using ([9]) and making the transformation z = — 7 ^^, 
we have 


poo poo 

E[Dn{g)] = Rniz,y)dzdy, 

Jo J —00 


Rn{z,y) = I(-^ \{z)Ki{z)fi{g{y) + zhi)f 2 {y). 


where 


From (USD and Assumptions 1 , 3 and 4, we have, for every hxed z and y, 
lim Rn{z,y) = Ki{z)fi{g{y))f2iy). 

n—^oo 

By Assumption 1, we have the dominance 

0 < Rn{z,y) < MfKi{z)f 2 {y). 


(15) 


(16) 


Assumption 1 and 2 ensure that the bounding function on the right hand side of 
(fTHp is integrable. Thus, by applying DCT we have 


lim 

n—>-oo 


poo 

E[Dn{g)] = / fi{g{y))f 2 {y)dy. 
Jo 


(17) 


Now from Lemma lA.ll proved below, we have lim^^oo Car [Il„(y)] = 0. This 
establishes m and completes the proof. □ 


Lemma A.l Under Assumptions 1, 2, 3 and 4, for any g G Go, we have 

lim Car[iV„(y)] = 0, 

n^oo 

lim Car[L»„(y)] = 0, 

n^oo 

where Nn{g) and Dn{g) are defined in (0) and (ED, respectively. 


(18) 

(19) 
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Proof of Lemma I A. 1 1 From (j8|). we have 
Var[Nn{9)) 


- nl n2 nl n2 ✓ / 

^ i=i j=iii=iji=i ^ 


-g(gj) 


{nin2hih2Y 


K2 


h2 


Ki 


ti' 

h. 


K2 


'Yt,-n, 

h2 


= P 1 + P 2 + F 3 + F 4 , 


( 20 ) 


where 


1 r //. 


V"2 = 


{nin2h 

1 

{nin2h 


-9{sj) 


2 = i J = 
nl n2 nl 


hi 


Yt. - E' 


1L\K2' 


(^0 


/12 ^ 

'Yu-Y^^ 

ho 




( 21 ) 


( 22 ) 




E 3 = 


(reinohi/io) 


'Ef. -E' ' 

St 


=i j=i/=i 

(^i) 

,,, 

nl n2 nl n2 


E 4 = 


{nin2hih2Y \ \ hi ) 

(^i) 


'Ef. -E' 

‘'l S-j 




ti'-g{sf)\ (Yt,, Y! 


hi 


K 2 


ho 


(23) 


We consider the convergence of each term on the right hand side of (|20p separately. 
By using ([T|) and (I^TI) . we have 


El = 


1 


Vii--^E‘^[Nn{g% 


where 

Ell = 


/ oo poo poo poo -I 

ir^ 

-ooJ—ooJo Jo ^1 


nin2hih2 nino 


1 t^ 2 f X- g{y) \ J_^2 (Mx) - migoiy)) + U-V 


hi } h 


xfi{x)f2{y)fei{u)fe2{v)dxdydudv. 
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By making the transformations zi = ^ and Z 2 = ™(go(^))+^ y ^ have 

/ OO POO POO POO 

/ / / -Jzi)K'^{zi)K^iz2)fi{giy) + zihi) 

-ooJo J-ooJ-oo ^ H 

xf 2 {y)fei{v - m{g{y) + zihi) + m{go{y)) + Z 2 h 2 )fe 2 {v)dzidz 2 dydv. (24) 

From Assumptions 1, 3, and 4 and the continuity of m, for any fixed real {zi,Z 2 ,y, v), 
a similar argument as given for m shows that the integrand function on the right 
hand side of (124)1 converges, as n —>• oo, to 

Kf{zi)K 2 {z 2 )fiigiy))f 2 {y)feii'^ - m{g{y)) + m{go{y)))f^ 2 {v). 


We have the dominance of the integrand function on the right hand side of 
by the integrable function 

M]Kf{zi)Kl{z2)f2{y)fe2{v). 

By using DCT and convergence of the integrand on the right hand side of (I24|) . we 
have 

/ OO POO 

Kl{zi)dzi / Kl{z2)dz2 

-OO j —OO 

/ OO POO 

/ fi{g{y))f 2 {y)feAv - m{g{y)) + m{go{y)))f^^{v)dydv. 

-OO J 0 


' —OO ^0 

Now, from (|14p . we have 


1 


nin2 


-E\Nn{g))=0 


Thus, we have 


Vi = 0 


1 


nin2hih2 


+ 0 


nin2 


nin2 


= O 


1 


nin2hih2 


We now consider the term V 2 . From (122)1 and ([T|), we have 


V"2 = 


ni-l 2 


nin2 


W 21 - 


nin2 


E\Nn{g)), 


(25) 


(26) 


Where 


^21 = 


X —K2 

n-2 

X ^Ki 
hi 


x-g{y) 

hi 


h2 

x' - g{y)\ 1 „ f m{x') - m{go{y)) + F - V 

vFf - h2 - 


/ OO POO POO POO POO POO 

-OO J—OO J—OO J 0 Jo Jo hi 
m{x) — rn{go{y)) + u — v 


X h{x)fi{x')f 2 {y)fei{u)fei{u')fe 2 {v)dxdx'dydudu'dv. 


(27) 
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By making the transformations zi = ^2 = ~ m(x) m(go(y))+^ v ^ 

and Zi = )-m{g^{y))+u -V ^ integrand on the right hand side of (1271) is 

I{-g{y)/hi,oo)izi)Ki{zi)I(^_g(^yyh^^oo)iz2)Ki{z2)K2{z3)K2{z4)fi{g{y) + Zihi) 

X fiiaiy) + Z 2 hi)f 2 {y)feAv - m{g{y) + zihi) + m{gQ{y)) + h 2 Z‘i) 

X fei{v - m{g{y) + Z2hi)+m{go{y))+ h 2 Zi)fe 2 {v). (28) 

From Assumptions 1, 3, and 4 and the continuity of m it follows via a similar 
argument for (fT 2 ]l that for any zi, 2 : 2 , 2 : 3 ,- 24 , y and v the above function converges, 
as n ^ 00 , to 

Ki{zi)Ki{z2)K2{z3)K2{zi)fl{g{y))f2{y)fl{v - m{g{y)) + m(yo(y)))/ 62 (^^)- 

The integrand function on the right hand side of (1281) is dominated by the integrable 
function 


A/j All ( 21 ) iFi ( 22 ) Ar 2 ( 2:3 ) A '2 (24 )/2 (y)/e 2 (^) ■ 


By using Assumption 2 and the convergence of 
right hand side of (l27)) . we have 


and applying DCT on the 


poo poo 

Urn V 21 = fU9iy))f2{y)fei{v-m{g{y))+m{go{y)))f^2{v)dvdy. 

Jo J -00 

Now, from (HH), the second term on the right hand side of ()26p turns out to be 


'JJ^E^Nnig)) = O 
nin2 


It follows that 


poo poo 

n2V2^ / / fi{g{y))f2{y)f^^{v-m{g{y)) + m{go{y)))fe2iv)dvdy-E^{Nn{g)) 

Jo J —00 


i.e., 


U = o(l). (29) 

By using a similar argument as for the term V 2 , we have 

poo foo f2( \ 

niV3^ / -^j^fi{g{y))feAv-m{g{y))+m{go{y)))fl{v)dvdy-E^{Nn{g)) 

Jo J -00 9 [y) 


i.e. 
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F3 = 0(i). (30) 

Finally, we consider the term V 4 . By using the model specification (HD, we have 

V 4 = 0. (31) 


The proof of (fTHjl is completed from (1251) . ([ 2 ^ . (I30]l . (f3Tl) and by using Assump¬ 
tions 3 and 4. 

We now compute Var[Dn{g)\- Note that, from ([9]), we have 
Var[Dn{g)\ = Ti + T 2 + T^ + T 4 , 


Where 
Ti = 

T 2 = 

n = 

n = 


nl n2 




^ j All 


ti-g{sj) 


^ nl n2 nl ^ y 


{nin2h 

1 


*=r i'=l 


nl n2 n2 

1 = 1 ] = l j' = l 

i^i) 


hi 

U -g{sj) 

hi 

tj - g{sj) 
/ii 


,Ki 


,Ki 


ti' gi^j) 

hi 

ti g{sj') 

h^ 


- nl n2 nl n2 y / 

^EEEEc- if. 

i=i j=ii'=ij'=i ^ 


{nin 2 h 


tj -g{sj) 


*=r 3=1 i'=lj'=l 


:Kl 


ti' gi^j') 


We consider the convergence of the terms Ti, T 2 , T 3 and T 4 separately. Consider 
the term Ti. By making the transformation 2 : = ^ and by using the model 
specification, we have 


where 


Ti 


1 

nin2hi 


Til 


1 

nin2 


E‘'[Dn{g)], 


(32) 


poo poo 

Tii= / I^_g(^yyhi,oo){z)Kf{z)fi{g{y) + hiz)f 2 {y)dzdy. (33) 

^0 J —00 


Note that the integrand on the right hand side of (l33|) is bounded by the in- 
tegrable function MfK‘l{z)f 2 {y)- Further, a similar argument, as used for the 
Var{Nn), shows that, for any given y and z, the integrand function converges to 
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Ki{z)fi{g{y))f 2 {y) as n —>■ oo. Thus, by applying DCT and Assnmptions 3 and 4, 
we have 

/ oo roo 

Kl{z)dz / fi{g{y))f 2 {y)dy. 

-oo Jo 

Now, from (dD, the second term on the right hand side of (1321) turns out to be 


Thus, we have 


nin2 


-E\Dn{g)]=0 


nin 2 


Ti =0 


1 


nin2hi 


+ 0 


nin 2 


= O 


1 


nin2hi 


We now consider the term T 2 . By making the transformations z = ^ and 
z' = we have from (fT|) 


T 2 


ni — 1 

nin 2 


T 21 


ni — 1 

nin 2 


E^[Dn{g)], 


(34) 


where 


roo roo 

Jo J- 


T 21 = 


g'(y)//ii,oo)(^)-^l(^)-^{—g(j/)//ii,oo)('2' )E\{^Z ) 


—00 J —00 


xfliaiy) + zhi)fi{g{y) + z'hi)f2{y)dzdz'dy. 


(35) 


Note that, the integrand on the right hand side of (1351) is bounded by the integrable 
function MjKi{z)Ki{z')f 2 {y)- A similar argument as nsed for the convergence 
of the term Tn shows that, for any given z,z' and y, the integrand function 
converges to Ki{z)Ki{z') {g{y)) f 2 {y), as n ^ 00 . Thus, by applying DCT and 
Assumption 2, 3 and 4, we have 

/•OO 

lim T 21 = / ff{giy))f 2 {y)dy. 

n-)-oo Jq 

Now, from (I17p . the second term on the right hand side of (|34p turns out to be 


ni — 1 

nin 2 


EHDnig)) 
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n2T2^ / fi{giy))f2{y)dy - E^{Dn{g)) 
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Thus, we have 
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i.e., 


T2=0 


1 

n2 


We now consider the term T^. A similar argument, as used for the convergence 
of the term T 2 , shows that 


niTs 


f 


9'{y) 


dy-E\Dn{g)) 


i.e., 


T 3 



Finally, by using the model specification, we have r 4 
proof of ([T9]). 


0. This completes the 
□ 


Proof of Theorem 13.31 Let us denote the convolution of the densities and 
/ej by fei+e 2 - We first show that fti+e 2 is strictly unimodal at zero. Indeed, for 
u > 0, we observe from Assumption 1 that 


/ CX) 

fei{v)fe2{u-v)dv 

-00 

POO POO 

= / fei{v)fe 2 {u-v)dv+ / /.j (u)/e 2 (u + u)d 

Jo Jo 


Because of the strict unimodality of 7^2, for U 2 > ui> 0, we have 


/ei+e2('^2) /ei+e2(’^l) 

POO 

= / feAv)[{fe2{u2-v) - fe2{ui-v)} + {fe^{u2+v) 

Jo 

< 0 . 


/ejlui +U)}] dv 


By a similar argument, the same inequality holds for U 2 < ui < 0. Thus, fei+e 2 
is strictly unimodal at zero. 

Now observe from ([T|) and Assumption 1 that 

, . _ ir fii9{y))f2{y)fedv - m{g{y)) + m{go{y)))f,^{v)dydv 

^ ir fii9iy))h{y)dy 

/o°° fii9iy))f2{y)fei+e2im{giy)) - m{go{y)))dy 
/o°° f 1(9(9)) f2{y)dy 
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In particular, L^qq) = Thus, 


L{9o) - L{g) 


f rni hi9{x))f2{x)fe^+e2{m{g{x))- 

’ h{g{x))h{x)dy 

/o°° h{g{x))f 2 {x) [Ai+caW - /,i+, 2 (m(g(x)) 
/o“ h{9{x))f2{x)dy 


m{go{x)))dx 


'iTx{ 9 o{x)))] dx 


Unimodality of /ei+ej 0 implies that 


[/€i+€ 2 ( 0 ) - fei+e 2 irn{gix)) - m{go{x)))] > 0 for all x. 


This inequality proves part (a). 

In order that the last expression for L{gQ) — L[g) happens to be zero for 
some S' G Go, the above difference must be equal to zero for all x such that 
/i(Sf(x))/ 2 (x) > 0, i.e., for x G S'g D Sg^, where Sg is as defined at the beginning of 
Section 3. Since fei+e 2 is strictly unimodal at 0, this requirement reduces to 


m{g{x)) = m{go{x)) 'ixeSgr\Sg^. 

It follows from Theorem 13.11 that Sg = Sg^ and g{x) = gQ{x) Vx G Sg^. This 
completes the proof of part (b). □ 


Proof of Theorem 13.41 In accordance with Assumptions 2 and 2A, let the 
positive real numbers c, Mk, be such that 0 < c < Arj(x) < Mk and |i7'(x)| < 
for f = 1,2. 

We first obtain a stochastic upper bound on the variation in L„(-). From d?]), 
for any given g,g G G, we have 


7'n(ff) Ln(^g)\ 


Nn{g) Nnjg) 

Dnig) Dn{g) 

\Nn{~g)Dn{g) - Nn{g)Dn{g)\ 

Dn{g)Dn{g) 


From ([9]) and Assumption 2A, we have Dn{g) > Therefore, 


l^nig) ^nig)] 

< ’^\Nn{9)Dn{g)-Nn{g)Dn{~g)\ 

& 

< ^ {B^{g) \Nn{g) - Nn{g)\ + Nn{~g) \Dn{~g) - Dni9)\} . 
& 


( 36 ) 
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We now compute the upper bounds for both the terms on the right hand side 
of Note that, from (l8|), we have 


\Nnig) - Nr^ig)\ 


^ nl n2 ^ 




nin2^ ^ hi 

1=1 j=i 
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rein2/ii/i2 


nl n2 
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h2 I 
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ti-g{sj)\} 1 (Yu -Ysj 


Ki 


i=l j=l 

By using the mean value theorem, we have 


hi 

ti-g{sj) 

hi 


-Ki 


h2 

tj - g{sj) 

hi 


\Nn{g) - Nn{g)\ 
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1 


nl n2 

EE^^-^ 


nin2hi/i2 ^ , 

1=1 j=i 

K[{xii{ti,Sj,g,g)) 


' y.-K \ 

g{sj) - g{sj) 


hi 


where 3:o(ti,Sj,5,5') ^ (^min 

from (|37l) and Assumption 2A, we have 


-9{sj) u-gjsi) 
hi ’ hi 


, max 


U-g{sj) U-gjsj) 

hi ’ hi 


\Nn{g) - Nn{g)\ < ^ • lb - 5ll • t^n, 


where 


Ur,. = 


1 


nl n2 


nin2h2 




ho 


(37) 
Now, 

(38) 

(39) 


i=i j=i 

We now turn to the second term on the right hand side of (I36|) . From Q, we have 

nl n2 


\Dn{g) - Dn{g)\ < 


1 


nin2hi 


EE 


Ki 


ti -9{sj)\ - 9{sj) 


h 


1 


-Ki 


i=l j=l 

From Assumption 2A and the mean value theorem, we have 

nl n2 


hi 


\Dn{g) - Dn{g)\ < 


1 


nin2hi 


EE K[{xi{ti,Sj,g,g)) 


i=i j=i 


9{sj) - g{sj) 


< 


M 




(40) 


, max 


wherewhere3:i(ti,Sj,5,5() G (^min 
by using (l38|l . (liil and (|36|) . we have 

\Ln{g) - Ln{g)\ < Bn{g)-\\g-g\\, 


hz|M,Ezg£i))y Now, 


(41) 
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where 

m' 

Bn{g) = {Nn{g) + Un ' I)n(ff)} • 

The expression on the right hand side of (I4ip gives an upper bound on the change 
of the functional Ln(-) with change in time transformation functions in G. 

Note that 


\Ln{g)-L{g)\ < \Ln{g)-Ln{g)\ + \L{~g)-Lig)\ + \Ln{~g)-L{~g)\, (42) 

where L(.) is defined as in Q- 

Set e > 0. Lemma Fa. 21 proved below, implies that there exists 5^ > 0 such that 


\\g - g\\ < implies \L{g) - L{g)\ < 


(43) 


Theorem l3.2l and Lemma lA.31 proved below, implies that for all g there exists Mg 

p 

such that Bn{g) Mg, which ensures 


P ( Bnig) > max <j —, 2Mg 


(44) 


Define Mg{g) = {g : \\g - g\\ < g} . For given g, let 


Kgp) = 


if M3>0 


if Mg = 0. 


(45) 


For g in Ms{^g^e){g)^ we have from (I^TI) 


\Lnig) - Lnig)\ < S{g,e) ■ Bn{g). 


(46) 


Note that {A/' 5 (g,,) ( 5 ) ■■ g € G} is an open cover of G. By Assumption 5, there 
exists a hnite sub-cover say \-^s{gj,e)igj)\ j with G C U^liA/'5(gj,e)(5j) for 

V J 7 = 1...Ale 

some finite k.. From 


7 ; 


' J = l...fee 

and (|36]l we have, 
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sup|L„( 5 ) - L{g)\ 
g€G 

< max sup iLn(g) - L(g)l 


< max 


sup \Ln{g) - Ln{gj)\ + sup \L{gj) - L{g)\ 


3-1,...,ke geMs(~g 


+ sup \Ln{gj) - L{gj)\ 

9G-'Vi(g )(gj) 


9&J^S(g.,e) iSj) 


J=^ 


{^(9j,e)Bn{gj) + I + \Ln{gj) - L{gj)\] 

€ 

max 6{gj,e)Bn{gj) + x + X] - L{gj)\- 

From (I45p and (I47p we have, 


i=i 


(47) 


P < sup|L„( 5 r) - L{g)\ > e 
Ugg 


< P 


^igjp)Bn{gj) > 31 + -P I \Ln{gj) - L{gj)\ > - 


Each summand of the first term on the right hand side of (1481) goes to zero by (1441) , 
while the second term goes to zero by Theorem 13.21 This completes the proof. □ 


Lemma A.2 Under Assumptions 1 and 5 the functional L{-) in & is uniformly 
continuous on G. 


Proof: Let 


/ OO /*oo 

/ fi{g{y))f 2 {y)fei{v-m{g{y))+ m{go{y)))fe 2 {v)dydv, (49) 

-OO Jo 
poo 

D{g)= / 

Jo 


fi{g{y))f2{y)dy. 


(50) 
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Then, L{g) = 

Let g € G and {gk G G; A: = 1, 2,...} be such that limfe_>.oo sup \gk — ^1 = 0. Now, 
from 

/ oo poo 

/ /i(s'fc(y))/2(y) 

-oo Jo 

X f^^{v - m{gk{y))+ m{go{y)))f^^{v)dydv. (51) 

Note that the integrand on the right hand side of (|5ip is bounded by the integrable 
function Mjf 2 {y)fe 2 iv)- Thus, applying DCT, we have 

/ oo poo f 

/ S lim /i(s'fe(y)) f/2(2/) 

-ooJo ) 

X lim f^^{v - m{gk{y)) + m(5o(y))) \ fe 2 {v)dydv. (52) 

P k^oo J 

Note that g^ ^ g ss k ^ oo pointwise. By using Assumption 1 and the fact that 
m is continuous, we have 


lim /^i(u 
k^oo 


lim fi{gk{y)) 

k^oo 

m{9k{y)) + rn{go{y))) 


fiiaiy)), 

feAv - m{g{y)) +m{go{y))). 


Thus, from (I52p . we have 


lim N{gk) 
k^oo 


/ oo poo 

/ fi{9{y))f2iy)feiiv 

-ooJo 

N{g). 


m{g{y)) + m{go{y)))fe 2 {v)dydv 


This shows that the functional N(-) is continuous on G. 

A similar argument shows that D{-) is also continuous. Further, note from 
Assumption 1 and (|5U|) that D{g) > 0 for any g G. This establishes that L is 
continuous on G. From Assumption 5, L is uniformly continuous on G. □ 


Lemma A. 3 Let 


^ nl n2 

nin2h2 ^ ^ 

1=1]=i 

Then, under Assumptions 1, 2, 3 and 4, 


If' — To. 

S'] 


( 53 ) 


/ OO poo poo 

/ / fi{x)f 2 {y)fei{v-m{x) + m{go{y)))f^^{v)dxdydv. 

-OO Jo Jo 
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Proof: From (j53p and ([T|), we have 


Em = r r r r 

J-oo J-oo Jo Jo ^2 V ^2 

X fi {x)f2 {y)fei {u)fe 2 {v)dxdydudv. 

By making the transformation w = ^ have 

/ oo /*oo /•oo /•oo 

/ / / K2iw)h{x)hiy) 

-OO J —oo ^0 ■/ 0 

X f^^{v - m{x) + m{go{y)) + wh 2 )fe 2 {v)dxdydwdv. 


(54) 


From Assumption 1, 3 and 4, the integrand on the right hand side of (I54p converges, 
as n —>■ oo, to 


K2{w)fi{x)f2{y)feAv - m(x) + rn{go{y)))f^2{v), 

and is bounded by the integrable function 

MfK2{w)fi{x)f2{y)fe2{v). 

By applying DCT, we have 

/ OO POO POO 

/ / fiix)f2{y) 

-oo Jo Jo 

X /e, (u - m(x) + m{go{y)))f^^ {v)dxdydv. (55) 

We now turn to the variance of C/„. From (|53p . we have 

Var{Un) = 4i + P 2 + + Kt, 


where 

Vh= 

F2= 

F3= 

P4=' 


ni 712 




'Yt- -YL 


ni 77,2 ni { 

\ ^2 


(nin2/i2)2^^tl^Xi I 

. ni 77,2 712 r 

I 

i^j) 

ni 712 Fil 712 


h2 Jj 

h2 

(Y.-Yl. 


,K2 


'Yt, -Y' 

S q 


h2 


1 


( ^2 


(¥^3) 


h2 ) 

( Yu-n^ 

I ho 


,K2 


'Yt^-YL 


:K2 


h2 y J ’ 

ho 


(56) 

(57) 


(58) 


• (59) 
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We show that the terms Vi, for i = 1,... ,4, converges to zero as n ^ oo. By 
making the transformation w = "^V)~"^( 9 o{y))+u-v (|56]1 and ([T]), we have 


where 




1 

nin2/i2 


^11 


1 

nin2 


E^Un), 


(60) 


/ OO poo POO POO 

/ / / Kliw)hix)f2iy) 

-OO J—oo 0 Jo 

X f^^{v - m{x) + m(5o(y)) + h 2 w)f^^{v)dxdydwdv. 


(61) 


From Assumption 1, 3 and 4, the integrand on the right hand side of (1611) converges, 
as n ^ oo, to 


K 2 {w)fi{x)f 2 {y)fei{v - m{x) + m{go{y)))fe 2 {v) 


and is dominated by the integrable function 

MfK^ {w)fi (x) /2 (y)/e 2 (v). 


Thus, by applying DOT, we have 

/ OO 

K2{w)dw 

-oo 


’ —oo 

/ oo POO POO 

/ / fi{x)f 2 iy)fei{v - m{x)+ m{goiy)))fe 2 iv)dxdydv. 

-OO Jo Jo 

Now, from (I55p . the second term on the right hand side of ()60p turns out to be 


1 


nin2 


-E\Un) = O 


1 


nin2 


Thus, we have 


Vi = 0 


1 


nin2/i2 


+ 0 


nin2 


= O 


1 


nin2h2 


By making the transformations w = ™(go(y))+^ v ^ _ m(x ) m(go(t /))+-» —v 

and using m and dSZI), we have 

1^2 = (62) 
nin2 nin2 

where 


/ OO POO POO POO POO POO 

/ / / / / K 2 iw)K 2 {w')hix)f 2 {y)fiix') 

-OO J—OO J—OO J0 Jo Jo 

X f^^{v-m{x) + m{go{y)) + h2w) 

X — m{x') + Tn{go{y)) + h 2 w')fe 2 iv)dxdydx'd'wdw'dv. (63) 
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From Assumption 1, 3 and 4, the integrand function on the right hand side of (I63p 
converges, as n —)• oo, to 


K2{w)K2{w')fi{x)f2{y)fi{x')fei{v - m{x) +m{go{y))) 
X f^^{v-m{x') + m{go{y)))fe 2 {v), 

and is dominated by the integrable function 


MjK2iw)K2{w')fi{x)f2{y)fi{x')fe2{v). 

Thus, by using DCT, we have 

/ OO /•OO ^ poo 2 

/ f2{y)fe2iv)\ / fi{x)f^^{v-m{x) + m{go{y)))dx\ dvdy. 

-ooJo Uo J 

Now, from ()55p . the second term on the right hand side of ()62h turns out to be 


ni — 1 

nin2 


E\Un) = 0{ — 


n2 


Thus, we have 

/ OO poo 

/ / 2 (y)/ 62 (^ 

-oo J 0 


i / h{x)feA'>^-'f^{x)+m{gQ{y)))dx\ dvdy 


/ OO poo poo A 2 

/ / h{x)f2{y)fei{v-m{x)+m{gQ{y)))fe2{v)dxdydv\ , 

-ooJO Jq j 


i.e., 


r. = 0| I, 


A similar argument, shows that 


/ OO poo f poo ^ 2 

/ fiix)fei{v)\ / f2iy)fe2{v + m{x)-m{go{y)))dy\ dvdx 
-OO ^0 1^0 J 

' poo poo poo 'j 2 

/ / / /i(a^)/2(y)/€i(t’-m(x)+m(5o(y)))/e2(^^)'^3:(iydu L 

. J —ooJ 0 0 J 


i.e.. 


^3 = 0 


1 

ni 


The term V 4 is seen to be 0 from the model specihcation. This completes the 
proof. □ 
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Proof of Theorem 13.5t For any given e > 0, we have 
P{\Ln{gn)-L{go)\>e} 

< P{\Ln{gn) - L{gQ)\ > e, |L„(fifo) - -^^(50)1 < e} + -P{|^n(5'o) - L{go)\ > e} 

^ PW^n{,gn) -^'( 50)1 ^ |-^'n(5o) -^(5o)| ^ |-^'n(5n) -^'(5n)| ^ c} 

+P{\Ln{gn) - L{gn)\ > e} + P{\Ln{go) - L{go)\ > e}, (64) 

where gn is as in ([3]) . We will complete the proof by establishing that all the three 
terms on the right hand side of (I64p are arbitrarily small. 

We begin with the first term on the right hand side of (j64|) . Note that, from Q, 
we have 

Pnigo) — Ln{gn)- 

Therefore, from Q, we have 

if \Ln{go) - L{go)\ < e then Ln(ffn) > L{go) - e. (65) 

We now turn to computing an upper bound for Ln{gn) in terms of L^go). From 
Theorem 13.31 we have L{gn) < L^go). Therefore, 

if \Ln{gn) - L{gn)\ < € then Ln{gn) < L{go) + e. (66) 

Further, and (fBUjl imply that 

if \Ln{go) - L{go)\ < € and |Tn(ffn) - Lign)\ < e then \Lnign) - L{go)\ < e. (67) 
Thus, from (fU71) . 

P{\Lnign) - L{9o)\ > e, {Lnigo) “ ^(ffo)| < e, \Ln{gn) “ L{gn)\ < e} = 0, 

which takes care of the first term on the right hand side of ()64p . 

We now consider the second term. Observe that 

\Ln{gn) - L{gn)\ < sup \Ln{g) - L{g)\. (68) 

g€G 

From (1681) and Theorem 13.41 we have 

p 

Ln{gn) P{,gn) ^0. 

This ensures that the second term on the right hand side of (|64p goes to zero as 
n —>■ oo. Further, Theorem (13.21) ensures that the last term on the right hand side 
of (f64)l goes to zero too. This completes the proof. □ 
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Proof of Theorem 13.6t If Qn ^ go, then there exists an e > 0 and a d > 0 such 
that 


P{sup \gn — go\ > e} > (5 infinitely often. (69) 

Note that, M^{go) = {g : sup 1^ — go\ > g ^ G} is a closed subset of G. 
From Assumption 5 and Lemma IA.21 there exists a g (z (g'o) such that g = 
argmaXgg_;\/-c(gp) L{g). It follows from part (b) of Theorem 13.31 that the supremum 
of the functional L is attained only at g^. Therefore, gn € Ff^igo) implies 

\L{go) - L{gn)\ = L{go) - L{gn) > L(c/o) - L{g) > 0. (70) 

Denote rj = L{gQ) — L{g). By using the triangular inequality, we have 

\L{go) - Ln{gn)\ + \Ln{gn) - L{gn)\ > \L{go) - L{gn)\. (71) 

From (fTOjl and (fTTT) . gn G M^{go) implies 

\L{go) - Ln{gn)\ > 

Now from (1721) . 

if gn G ^feigo) then sup \Lnig) - L{g)\ 
geG 

Therefore, from (I73p . we have 

p!^\Ln{gn)-L{go)\ > || 

> P< G M^{go) and sup 

I geG 

> P{sup\gn - go\ > e} + P 

From (16911 . the first term on the right hai 
often. From Theorem l3.4l the second term on the right hand side of (I74p is greater 
than 1 — I for all but finitely many n. Therefore, 

P {\Lnign) - L{go)\ >1} > — infinitely often. 


T] - \Ln{gn) - L{gn)\- (72) 

< I implies \Lnign) - L{go)\ > (73) 

|L„(9)-i(9)|<|}, 

I. sup\Lnig) - L{g)\ - 1. (74) 

IseG 4 J 

nd side of (I74p is greater than 6 infinitely 


This contradicts Theorem 13.51 and completes the proof. 


□ 
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Table 1: Average normalized IMSE (in 10 


Method 

Seen. 1 

Seen. 2 

Seen. 3 

Seen. 4 

Cont. mon. registration 

1.000 

1.022 

0.598 

0.108 

Self-modelling registration 

0.375 

- 

0.433 

- 

The proposed method 

0.104 

0.287 

0.179 

0.189 


Table 2: Some descriptive statistics 


Data 

Data set 1: Vostok 

Data set 2: EPICA Dome 

Size 

Range (T -Value) 

Size 

Range (V -Value) 

Carbon dioxide 

283 

182.2-298.7 

537 

183.8-298.6 

Methane 

457 

318-773 

1545 

342-907 

Temp, deviations 

3,310 

(-)9.39-3.23 

5028 

(-)10.58-5.46 


Table 3: Average squared difference between pair of data sets 



Carbon dioxide 

Methane 

Temp. dev. 

Pre-alignment 

266.09 

3739.06 

3.04 

Post-alignment: 




Cont. mon. registration 

70.76 

- 

1.76 

Proposed method 

16.20 

1353.93 

1.31 













Carbon dioxide Concentration (parts per million voiume) 
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Figure 1: Ice core data on the atmospheric concentration of carbon dioxide 
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Figure 2: Point-wise bias, standard deviation, and MSE of the estimators 
of Qq by continuous monotone registration (broken line), self-modelling reg¬ 
istration (dotdash) and the proposed method (solid) under scenario 1 (top 
row), scenario 2 (second), scenario 3 (third), and scenario 4 (bottom) 
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Figure 3: Alignment of data sets on atmospheric concentration of carbon 
dioxide 










Temp, deviations Temp, deviations Temp, deviations 
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Figure 4: Alignment of data sets on average temperature deviations 






Methane concentration Methane concentration 
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Figure 5: Alignment of data sets on atmospheric concentration of methane 
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SUPPLEMENTAL MATERIALS 


(i) gO as in (5) 


(ii) gO as in (6) 




Figure SI: Qq functions for simulations 
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(iii) Methane 
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Figure S2: Estimates of go for paleoclimatic data sets 







